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A  finite-field  method  for  the  calculation  of  polarizabilities  and  hyperpolarizabilities  is  developed  based 
on  both  an  energy  expansion  and  a  dipole  moment  expansion.  This  procedure  is  implemented  in  the 
MOPAC  semiempirical  program.  Values  and  components  of  the  dipole  moment  (p),  polarizability  (a), 
first  hyperpolarizability  and  second  hyperpolarizability  fy)  are  calculated  as  an  extension  of  the 
usual  MOPAC  run.  Applications  to  benzene  and  substituted  benzenes  are  shown  as  test  cases  utilizing 
both  MNDO  and  AMI  Hamiltonians. 


INTRODUCTION 

The  polarization,  P,  induced  in  a  medium 
by  an  external  electric  field  F  is  given  by 

P  =  P0  +  X‘"-F  +  x'2'FF 

+  x'3,  •  F  ■  F  ■  F  +  ...  (1) 

where  y'"'  is  the  /?th  order  susceptibility  ten¬ 
sor  of  the  bulk  medium.1  The  bulk  suscepti¬ 
bilities  can  in  turn  be  expressed  in  terms  of 
the  molecular  induced  dipole.  The  dipole  mo¬ 
ment  of  a  molecule  interacting  with  an  elec¬ 
tric  field  can  be  written 

Mi  ~  Mi’  +  otjjFj  +  (1/2 )(i,jkFjFk 

+  (l/6)y  ijkiFjFkFi  +  ...  (2) 

where  p"  is  the  permanent  dipole  moment 
and  a,r  (3ijk.  and  y,ikl  are  tensor  elements  of 
the  linear  polarizability  and  the  first  and 
second  hyperpolarizabilities,  respectively,  of 
the  molecule.1  The  second-order  term  gives 
rise  to  sum  and  difference  frequency  mixing 
(including  second  harmonic  generation)  and 
optical  rectification.  The  third-order  term  is 
responsible  for  third  harmonic  generation  and 
two-photon  resonances.  Typical  applications 
where  these  nonlinear  effects  are  important 
include  amplification,  high-resolution  spec¬ 
troscopy,  picosecond  pulse  generation,  infra¬ 
red  image  conversion,  image  transmission 
through  optical  fibers,  optical  processing,  and 
information  transfer  in  general.2,3 

As  part  of  a  long-range  project  to  calculate 
the  nonlinear  optical  properties  of  polymer 


and  large  organic  systems,  procedures  for  the 
calculation  of  molecular  hyperpolarizabilities 
have  been  implemented  into  the  well  estab¬ 
lished  semiempirical  electronic  structure 
program  MOPAC.4  In  this  article,  the  details 
of  the  method  of  calculation  will  first  be  out¬ 
lined  and  then  the  method  evaluated  by  study¬ 
ing  a  series  of  substituted  benzenes. 


THEORY  FOR  CALCULATING 
POLARIZABILITIES: 

The  energy  of  a  system  in  an  electric  field 
can  be  written  as 

F(F)  -  £7(0)  -  p,F,  -  (1/2! )aijFiFJ 
- 

-  n/4'.)yljklF, FjFkF,  -  ....  (3) 

The  above  equation  contains  implied  sums 
over  repeated  indices,  E( 0)  is  the  energy 
with  no  field  present,  and  F,  are  the  compo¬ 
nents  of  the  applied  field. 

If  the  molecule  is  considered  to  be  in  a 
uniform  electric  field  aligned  along  one  of 
the  axis  of  the  system  (i.e.,  IF,,  0,01),  the 
values  of  the  polarizabilities  along  that  axis 
<ju,, and  yxxxx )  can  be  obtained.  For 
this  case  the  energy  expression  reduces  to 

E(F, )  =  £7(0)  -  p,F,  -  (1/2 )auFf 

.  -  (l/6)A(1Ff  -  <l/24)y„(IF4  -  .... 

(4) 
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Truncating  this  expression  after  the  F4  term 
and  evaluating  the  energy  at  four  field 
strengths  (±F, ,  ±2F, )  leads  to  four  equa¬ 
tions  in  four  unknowns.  These  can  be  solved 
for  explicitly  and  the  results  are 

/x,F,  =  -(2/3) [F(F, )  -  F(-F,)J 

+  (1/12) [F(2F, )  -  E(~2F,)]  (5a) 

ot,,Ff  =  (5/2)F(0)  -  (4/3)  [F(F,)  +  E(-F,)] 

+  (1/12)  [F(2F, )  +  F(  -2F, )]  (5b) 

A„F?  =  [E(F  i )  -  F(-F,)] 

-  (1/2)[F(2F,)  -  F(-2F,)]  (5c) 

It  is  also  important  to  obtain  the  "nonaxial” 
components  of  the  polarizabilities.  For  ex¬ 
ample,  a,j  is  necessary  to  the  verify  rota¬ 
tional  invariance  of  the  results  and  to  Find 
the  principal  "optical”  axes.  The  values  of 
(iw  and  y,w  are  needed  for  a  comparison  with 
experimental  quantities.  To  obtain  the  com¬ 
ponents  with  only  two  different  indices  (i 
and  j),  a  linear  electric  field  is  applied  in  the 
ij  plane  along  45°  lines  between  axis  i  and 
axis  j.  In  this  case  the  energy  expression  is 

E(F,,Fj)  =  E( 0)  -  ix, F,  -  iijFj  -  (1/2 )auF‘f 

-  (1/2  )otjjFj  -  ctjjF'Fj 

-  (1  /6)/3,„F?  -  (1  /6)pjjjF] 

-  (l/2)pijjFlFj  -  (1  /WjuFjF? 

-  (l/24)y ,,„F4  -  (1/24)7 'mFj 

-  (1/6)7 ijjjFjF]  -  (1/6)7 auFjFf 

-  (1/4)7 iWFf  (6) 

If  the  "axial”  calculations  have  already 
been  done,  then  there  are  only  six  unknowns 
in  this  expression.  By  using  the  results  of 
energy  calculations  at  (F, ,  F,),  (F, ,  -F,), 
(~F,,Fj),  ( —  F, ,  -Fj),  it  is  possible  to  solve 
for  four  of  the  unknowns.  The  results  of  the 
manipulations  of  the  energy  expressions  are 

a,jF,Fj  =  (1/48)[F(2F„2F,)  -  E(2Fi,-2FJ) 
-  E(-2F,,2Fj) 

+  E(-2F„  -2 Fj)] 

-  (1/3) [E{F,,Fj)  -  E{F, ,  -Fj ) 

-  E(-F,,F}) 

+  E(  — F, ,  —Fj )  (7  a) 

PVJF,Fj  =  (1/2 )[E(-Ft,-Fj)  -  E(Fj,Fj) 

+  E(-F,,Fj)  -  E(Fi ,  -Fj )] 
+  IF(F,)  -  F(-F,)J  (7b) 


7 ujjF'fF'f  -  -4F(0)  -  [F(F , ,F; ) 

+  Ei  —  F, ,  —  Fj ) 

+  F(F , ,  -F , ) 

+  E(  -F ,F, )  J 

+  2[F(F, )  +  E(  -F, )] 

+  2[F(F;)  +  F(-F;)J  (7c) 

Further  calculations  must  be  done  to  solve 
for  7yy  and  yjm.  The  "axial”  equations  are 
the  same  as  originally  derived  by  Bartlett 
and  Purvis.3 

An  alternate  method  for  obtaining  the  po¬ 
larizability  and  hyperpolarizabilities  is  to 
use  the  equation  for  the  induced  dipole  mo¬ 
ment  (eq.  (2))  instead  of  the  energy.  By  using 
a  similar  procedure  for  evaluating  the  dipole 
moment  at  various  field  strengths,  the  fol¬ 
lowing  equations  can  be  derived 

ix,  =  (2/3)  [p.,(F,)  +  ix,(-F,)] 

-  (1/6) |ju.,(2F, )  +  ix,(-2F,)] 

(8a) 

anF,  =  (2/3)  [/x, (F, )  -  ix, (-F, )] 

-  ( 1/12)  [/x,  (2F, )  -  ix,(-2F, )] 

(8b) 

a.jFj  =  (2/3) [/*, (F,)  -  fXji-Fj)] 

-  (1/12)Lu,(2F/)  -  nS~2 Fj)] 

(8c) 

P.nF'f  =  (1/3)[m,(2F,)  +  nd-2  F.) 

-  fx,(F, )  -  iXj(-Fj)]  (8d) 

PwFj  =  (1/3)  [/a,  (2F,)  +  Vi(-2Fj) 

-  fXi(Fj)  -  ix, ( —Fj )]  (8e) 

7„„Ff  =  (1/2)  [fi,(2F , )  -  n,{-2F,)] 

~  lv-i(F,)  -  ix,(-F ,)]  (8f) 

7 i,jjFtF'j  =  (1/2 )[fx,{Fi,Fl)  -  ix, ( -F , , F, ) 

+  ix, (F, ,  -Fj) 

~  Hi(-F,,-Fj)] 

-  [fx,  (F , )  -  ix,  (-F,)]  (8g) 

Although  expressed  quite  differently,  these 
expressions  are  equivalent  to  those  recently 
published  by  Williams.6  The  main  differ¬ 
ences  between  this  work  and  Williams’  are 
in  the  number  of  calculations  that  must  be 
done  and  the  electric  fields  used. 

The  numerical  accuracy  of  the  finite  field 
equations  developed  is  sensitive  to  the  preci¬ 
sion  in  the  energy  or  dipole  moment  calcu¬ 
lation.  For  example,  it  can  be  seen  from 
eq.  (5d)  that  for  a  field  strength  of  0.001  au, 
an  error  of  about  10" 15  in  the  energies  will 
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lead  to  a  10"3  error  in  y.  For  a  similar  error 
in  the  dipole  moment  as  smaller  error  will 
result  in  y  because  of  the  lower  power  of  the 
field.  If  however,  in  trying  to  reduce  this 
error  by  using  too  large  a  value  of  the  field, 
then  electron  configurational  changes  may 
take  place  and/or  the  higher  terms  in  the  fi¬ 
nite  field  equations  will  become  important. 
In  a  subsequent  section  a  numerical  study 
will  be  given  to  demonstrate  the  importance 
of  these  effects. 

Both  methods  (the  energy  expansion  and 
the  dipole  expansion)  have  been  implemented 
in  the  MOPAC  program.  Since  the  results 
based  on  SCF  calculations  should  be  identi¬ 
cal  for  both  methods,  this  provides  a  good 
check  on  the  results  as  to  when  numerical 
errors  or  configuration  changes  take  place. 
If  Cl  calculations  are  performed,  the  results 
now  differ.  This  is  because  an  arbitrary  Cl 
wave  function  does  not  satisfy  the  Hellmann- 
Feynman  Theorem. 

Once  all  the  components  of  the  hyperpo¬ 
larizabilities  have  been  obtained  it  is  neces¬ 
sary  to  convert  them  into  "experimental” 
quantities.2  3  For  the  second  hyperpolariz¬ 
ability,  the  quantity  of  interest  is  the  mean 
value  given  by 

y  =  l/5{yxxxx  +  Xvvvy  +  Yzzzz 

+  2[y«vy  +  Yxxzz  +  Xvy«]} 

For  the  first  hyperpolarizability,  /3,  the  vec¬ 
tor  component  along  the  dipole  moment  is 
important.  The  quantity  of  interest  is 

a,  =  (3/5){0-/*}/iy  do) 

where  (3  ■  //.  =  (3xfix  +  pyfiy(3  +  (3z/iz  and 

Pi  =  {Pm  +  Pw  +  Pm }  (ID 

BENZENE 

Initial  calculations  were  carried  ouc  ^n 
benzene,  both  as  a  base  for  studying  the  sub¬ 
stituted  benzenes  and  as  a  well-studied  bench¬ 
mark.  Since  this  molecule  has  a  center  of 
symmetry,  the  dipole  (fx)  and  first  hyperpo- 
larizabilitv  ((3)  are  zero.  MNDO  based  values 
obtained  for  the  other  constants  (a  and  y) 
are  shown  in  Table  I.  It  should  be  noted  that 
in  the  calculation  of  the  polarizability  (a), 
the  atomic  correction  factors  developed  by 
Dewar  and  Stewart7  are  used.  The  use  of 
these  factors  gives  extremely  good  agreement 
with  experimental  results  for  a.  The  value 


Table  I.  Benzene  reults. 


Polarizability 

Second  hyperpolarizability 

XX 

89.06 

XXXX 

1835.7 

vv 

89.06 

vvvv 

1835.7 

zz 

28.17 

zzzz 

16.7 

xxw 

612.0 

«<au) 

68.43 

xxzz 

532.4 

at  A1) 

10.14 

yyzz 

532.3 

ytau) 

1408.3 

y(  10  <h  esui 

0.71 

Exp." 

10.33 

Exp." 

0.98 

'Dewar  and  Stewart,  1984. 
b  Reference  8. 


for  y  is  in  good  agreement  with  the  result  of 
Williams6  and  the  recent  experimental  esti¬ 
mate  of  the  static  field  value.8 

It  should  be  pointed  out  that  to  obtain  the 
type  of  numerical  precision  shown  in  Table  I 
( y xxxx  ~  Xvvvy  nnd  y xxzz  —  y yyZZ )  a  very  small 
SCF  convergence  criterion  must  be  used.  In 
this  work  a  value  of  1.0  x  10  20  was  typi¬ 
cally  used  —  much  smaller  than  is  normal 
for  semiempirical  energy  and  structure 
calculations. 

MONO-SUBSTITUTED  BENZENES 

In  order  to  further  evaluate  the  computa¬ 
tion  of  hyperpolarizabilities,  a  series  of  mono- 
substituted  benzenes  were  studied.  As  these 
molecules  do  not  possess  a  center  of  symme¬ 
try,  they  have  nonzero  values  for  the  dipole 
moment  and  are  good  test  cases  for  the  first 
hyperpolarizability,  (3.  A  comparison  of  the 
values  obtained  and  experimental  values 
are  given  in  Table  II.  There  does  not  seem 
to  be  any  systematic  differences  between 
theory  and  experiment,  particularly  when 
the  large  variation  in  experimental  results 
is  considered. 

The  results  are  seen  to  be  quite  good  for  F, 
CN,  and  OH,  with  essentially  the  same  re¬ 
sults  obtained  with  the  MNDO  and  AMI 
methods.  In  each  of  these  calculations  the 
optimized  structure  for  each  method  was 
employed.  However,  for  aniline  and  nitro¬ 
benzene,  the  situation  is  different.  For  aniline, 
the  AMI  results  agree  well  with  experiment 
but  the  MNDO  results  are  too  small  and  of  a 
different  sign.  However,  if  the  calculations 
are  carried  out  with  the  MNDO  method  at 
the  AMI  geometry,  a  value  of  0.66  x  10  30  esu 
is  obtained;  but  if  AMI  is  used  at  the  MNDO 
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Table  II. 


F 

CN 

OH 

NO, 

NH, 

Dipole  moment  (/x,  Debye) 
MNDO 

1.956 

3.342 

1.164 

5.187 

1.464 

AMI 

1.573 

3.336 

1  233 

5.240 

1.542 

Exp.  (Reference  9) 

1.66 

3.93 

1.55 

4.4 

1.53 

Hvperpolarizability  if}„,  10” 
MNDO 

’  esu) 

-0.65 

0.43 

-0.40 

0.39 

-0.19 

AMI 

-0.70 

0.41 

-0.51 

0.02 

0.83 

Exp.a 

0.44 

0.48 

0.17 

2.0 

0.79 

(Reference  10) 

0.53 

0.36 

2.3 

0.89 

0.70 

2.2 

1.23 

1.48 

11  Experimental  values  are  absolute  values  of  /3M. 


geometry,  the  result  is  0.01  x  10  30  esu.  This 
clearly  indicates  that  the  differences  are 
mainly  due  to  the  optimized  geometry  from 
each  method.  T  ie  better  AMI  results  agree 
with  the  experience  that  AMI  geometries  are 
better  than  MNDO  geometries. 

For  nitrobenzene  neither  of  the  methods 
give  correct  results.  For  this  molecule  the 
calculations  were  both  done  at  the  AMI  ge¬ 
ometry  because  of  an  unrealistic  orientation 
of  the  N02  group  perpendicular  to  the  ring 
in  the  MNDO  optimized  geometry.  The  prob¬ 
lems  with  nitrobenzene  do  not  seem  to  be 
problems  with  the  finite  field  methods  but 
with  the  parameterization  of  N  in  the  MNDO 
and  AMI  approximations.  This  may  also  ac¬ 
count  for  some  of  the  discrepancy  found  in 
aniline.  A  new  version  of  MOPAC  contain¬ 
ing  reparameterized  atoms  has  been  devel¬ 
oped  and  will  be  applied  to  this  problem.11 


The  question  of  the  numerical  stability  and 
the  "best”  field  strengths  to  use  in  eqs.  (5), 
(7),  and  (8)  is  addressed  by  the  results  shown 
in  Figures  1  and  2.  These  figures  show 
the  variation  of  the  calculated  values  of  (3 
and  y  for  fluorobenzene  with  field  strength 
(from  .0001  to  .1).  Over  a  wide  range  of  field 
strength,  the  energy  and  dipole  expressions 
can  be  seen  to  give  the  same  results.  For 
y  (Fig.  2),  erratic  values  are  found  below 
0.001  au  and  are  due  to  the  onset  of  numeri¬ 
cal  errors.  At  higher  field  strengths  the  com¬ 
puted  value  of  y  rises  until  a  dramatic  change 
occurs  between  0.1  and  0.05  au.  The  slow 
increase  in  magnitude  is  also  found  for  (3 
and  can  be  attributed  to  the  importance  of 
higher  order  terms  that  were  neglected  in 
the  expansion. 

To  further  evaluate  the  results  for  the  AMI 
and  MNDO  approximations,  a  detailed  com- 
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Figure  1.  Variation  of  /3  with  ‘ield  strength  for  fluorobenzene.  Filled 
circles  were  obtained  using  th  „■  energy  based  equations  and  open  circles 
were  obtain  with  the  dipole  based  equations. 
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circles  were  obtained  using  the  energy  based  equations  and  open  circles 
were  obtain  with  the  dipole  based  equations. 


parison  of  individual  components  calculated 
for  /j.  and  (i  is  shown  in  Table  III.  The  mole¬ 
cules  are  all  oriented  such  that  the  principal 
moments  of  inertia  lie  along  the  x,  y,  and  z 


axes.  This  comparison  reveals  that  for  ani¬ 
line,  the  agreement  is  quite  good  between 
the  two  methods  with  the  exception  of  the  fix 
and  /3XXX. 


Table  III.  Calculated  Components  of  /a  and  /3 


F 

CN 

OH 

NO, 

NH, 

Dipole  Moment  l/ai 
MNDO 

X 

0.770 

-1.315 

0.102 

-2.041 

-0.062 

V 

0.0 

0.0 

0.447 

0.0 

0.002 

z 

0.0 

0.0 

-0.001 

0.00 

0.573 

AMI 

X 

0.619 

1.313 

0.179 

-2.062 

-0.306 

V 

0.0 

0.0 

0.451 

0.0 

0.0 

z 

0.0 

0.0 

0.0 

0.0 

0.524 

Hyperpolarizability  (£0 

MNDO 

xxx  -135.19 

116.59 

-274.21 

-165.81 

-378.67 

JCVV 

8.41 

-22.19 

21.74 

90.37 

23.69 

xzz 

1.77 

-12.28 

-1.36 

1.19 

-5.73 

vvv 

0.0 

0.0 

-5.88 

0.0 

-0.14 

yxx 

0.0 

0.0 

-12.23 

0.0 

-0.24 

yzz 

0.0 

0.0 

-2.43 

0.0 

0.05 

zzz 

0.0 

0.0 

0.00 

0.0 

-2.39 

zxx 

0.0 

0.0 

-0.03 

0.0 

-66.26 

zyy 

0.0 

0.0 

0.01 

0.0 

-6.46 

AMI 

xxx  - 

144.12 

112.57 

-231.81 

-81.40 

-455.69 

xvv 

6.71 

-21.39 

20.15 

78.48 

28.43 

xzz 

1.86 

-12.89 

-0.71 

-0.16 

-8.87 

vvv 

0.0 

0.0 

-8.00 

0.0 

0.03 

yxx 

0.0 

0.0 

-10.83 

0.0 

0.04 

yzz 

0.0 

0.0 

-2.44 

0.0 

0.00 

zzz 

0.0 

0.0 

0.00 

0.0 

-4.54 

zxx 

0.0 

0.0 

0.00 

0.0 

-58.62 

zyy 

0.0 

0.0 

0.00 

0.0 

-5.36 
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CONCLUSIONS 

The  main  goal  and  strength  of  the  proce¬ 
dure  outlined  here  is  that  everything  is  done 
by  one  program:  geometry  optimization  and 
polarizabilities.  All  that  is  needed  is  a  rough 
initial  structure.  Since  the  development  is 
based  on  semiempirical  methods,  treatment 
of  large  organic  molecules  and  polymers  is 
possible.  While  large  basis  set  ab  initio  cal¬ 
culations  can  undoubtedly  do  a  more  accu¬ 
rate  job  on  small  molecules,  these  methods 
are  not  now  feasible  for  routine  calculations 
that  are  necessary  for  screening  molecules 
for  nonlinear  activity. 

One  important  extension  of  this  work  is 
to  test  the  usefulness  of  MNDO  and  AMI 
methods  to  calculate  frequency  dependent 
hyperpolarizabilities.  Work  is  currently  under¬ 
way  to  implement  this  via  a  sum-over-states 
approach  which  as  been  shown  to  work  quite 
well  within  a  CNDO  framework.12  The  fre¬ 
quency  dependent  values  are  necessary  to 
directly  compare  with  the  experimental  val¬ 
ues  and  to  make  judgements  of  the  potential 
usefulness  of  new  products. 

Research  sponsored  by  the  Air  Force  Office  of  Scien¬ 
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Contract  F49620-85-0-0013. 
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